home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / SOLUTION.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  4KB  |  144 lines
  1. //$$ solution.cpp                    // solve routines
  2.  
  3. // Copyright (C) 1994: R B Davies
  4.  
  5.  
  6. #define WANT_STREAM                  // include.h will get stream fns
  7. #define WANT_MATH                    // include.h will get math fns
  8.  
  9. #include "include.h"
  10. #include "boolean.h"
  11. #include "myexcept.h"
  12.  
  13. #include "solution.h"
  14.  
  15. int SolutionException::action = 1;
  16.  
  17. int iabs(int i) { return (i>=0) ? i : -i; }
  18.  
  19. SolutionException::SolutionException(char* c) : Exception(iabs(action))
  20. {
  21.    if (action) cout << "Error detected by solution package\n"
  22.       << c << "\n";
  23.    if (action < 0) exit(1);
  24. };
  25.  
  26. inline Real square(Real x) { return x*x; }
  27.  
  28. void OneDimSolve::LookAt(int V)
  29. {
  30.    lim--;
  31.    if (!lim) Throw(SolutionException("Does not converge"));
  32.    Last = V;
  33.    Real yy = function(x[V]) - YY;
  34.    Finish = (fabs(yy) < acc);
  35.    y[V] = vpol*yy;
  36. }
  37.  
  38. void OneDimSolve::HFlip() { hpol=-hpol; State(U,C,L); }
  39.  
  40. void OneDimSolve::VFlip()
  41.    { vpol = -vpol; y[0] = -y[0]; y[1] = -y[1]; y[2] = -y[2]; }
  42.  
  43. void OneDimSolve::Flip()
  44. {
  45.    hpol=-hpol; vpol=-vpol; State(U,C,L);
  46.    y[0] = -y[0]; y[1] = -y[1]; y[2] = -y[2];
  47. }
  48.  
  49. void OneDimSolve::Linear(int I, int J, int K)
  50. {
  51.    x[J] = (x[I]*y[K] - x[K]*y[I])/(y[K] - y[I]);
  52.    // cout << "Linear\n";
  53. }
  54.  
  55. void OneDimSolve::Quadratic(int I, int J, int K)
  56. {
  57.    // result to overwrite I
  58.    Real YJK, YIK, YIJ, XKI, XKJ;
  59.    YJK = y[J] - y[K]; YIK = y[I] - y[K]; YIJ = y[I] - y[J];
  60.    XKI = (x[K] - x[I]);
  61.    XKJ = (x[K]*y[J] - x[J]*y[K])/YJK;
  62.    if ( square(YJK/YIK)>(x[K] - x[J])/XKI ||
  63.       square(YIJ/YIK)>(x[J] - x[I])/XKI )
  64.    {
  65.       x[I] = XKJ;
  66.       // cout << "Quadratic - exceptional\n";
  67.    }
  68.    else
  69.    {
  70.       XKI = (x[K]*y[I] - x[I]*y[K])/YIK;
  71.       x[I] = (XKJ*y[I] - XKI*y[J])/YIJ;
  72.       // cout << "Quadratic - normal\n";
  73.    }
  74. }
  75.  
  76. Real OneDimSolve::Solve(Real Y, Real X, Real Dev, int Lim)
  77. {
  78.    Tracer et("OneDimSolve::Solve");
  79.    lim=Lim;
  80.    if (Dev==0.0)Throw(SolutionException("Dev is zero"));
  81.    L=0; C=1; U=2; vpol=1; hpol=1; y[C]=0.0; y[U]=0.0;
  82.    if (Dev<0.0) { hpol=-1; Dev = -Dev; }
  83.    YY=Y;                                // target value
  84.    x[L] = X;                            // initial trial value
  85.    LookAt(L); if (Finish) goto finish;
  86.    if (y[L]>0.0) VFlip();               // so Y[L] < 0
  87.    x[U] = X + Dev * hpol;
  88.    LookAt(U); if (Finish) goto finish;
  89.    if (y[U] > 0.0) goto captured1;
  90.    if (y[U] == y[L])
  91.       Throw(SolutionException("Function is flat"));
  92.    if (y[U] < y[L]) HFlip();             // Change direction
  93.    State(L,U,C);
  94.    for (i=0; i<20; i++)
  95.    {
  96.       // cout << "Searching for crossing point\n";
  97.       // Have L C then crossing point, Y[L]<Y[C]<0
  98.       x[U] = x[C] + Dev*hpol;
  99.       LookAt(U); if (Finish) goto finish;
  100.       if (y[U] > 0) goto captured2;
  101.       if (y[U] < y[C])
  102.          Throw(SolutionException("Function is not monotone"));
  103.       Dev *= 2.0;
  104.       State(C,U,L);
  105.    }
  106.    Throw(SolutionException("Can't locate a crossing point"));
  107.  
  108. captured1:
  109.    // cout << "Captured - 1\n";
  110.    // We have 2 points L and U with crossing between them
  111.    Linear(L,C,U);                   // linear interpolation - result to C
  112.    LookAt(C); if (Finish) goto finish;
  113.    if (y[C] > 0.0) Flip();            // Want y[C] < 0
  114.    if (y[C] < 0.5*y[L]) { State(C,L,U); goto binary; }
  115.  
  116. captured2:
  117.    // cout << "Captured - 2\n";
  118.    // We have L,C before crossing, U after crossing
  119.    Quadratic(L,C,U);                // quad interpolation - result to L
  120.    LookAt(L); if (Finish) goto finish;
  121.    if ((x[L] - x[C])*hpol <= 0.0 || (x[L] - x[U])*hpol >= 0.0)
  122.       { State(C,L,U); goto captured1; }
  123.    State(C,L,U);
  124.    // cout << "Through first stage\n";
  125.    if (y[C] > 0.0) Flip();
  126.    if (y[C] > 0.5*y[L]) goto captured2;
  127.    else { State(C,L,U); goto captured1; }
  128.  
  129. binary:
  130.    // We have L, U around crossing - do binary search
  131.    // cout << "Binary\n";
  132.    for (i=3; i; i--)
  133.    {
  134.       x[C] = 0.5*(x[L]+x[U]);
  135.       LookAt(C); if (Finish) goto finish;
  136.       if (y[C]>0.0) State(L,U,C); else State(C,L,U);
  137.    }
  138.    goto captured1;
  139.  
  140. finish:
  141.    return x[Last];
  142. }
  143.  
  144.